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Abstract. I use the method of classical density- functional theory in the weighted- 
density approximation of Tarazona to investigate the phase diagram and the interface 
structure of a two-dimensional lattice-gas model with three phases - vapour, liquid, 
and triangular solid. While a straightforward mean-field treatment of the interparticle 
attraction is unable to give a stable liquid phase, the correct phase diagram is obtained 
when including a suitably chosen square-gradient term in the system grand potential. 
Taken this theory for granted, I further examine the structure of the solid-vapour 
interface as the triple point is approached from low temperature. Surprisingly, a 
novel phase (rather than the liquid) is found to grow at the interface, exhibiting 
an unusually long modulation along the interface normal. The conventional surface- 
melting behaviour is recovered only by artificially restricting the symmetries being 
available to the density field. 
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1. Introduction 

The statistical behaviour of a spatially inhomogeneous classical fluid (as is the case of the 
interface between two coexisting bulk phases or that of a fluid in a confined geometry) is 
best analysed within the (classical) density-functional theory (DFT) [HE]. Nowadays, 
the DFT method is widely recognized as being one of the most accurate theoretical tools 
in statistical mechanics, as also proved by the increasing number of applications in the 
scientific literature |3]. 

Generally speaking, the DFT provides a functional relation between the 
thermodynamic grand potential of the system and its local-density profile n(x). In 
practice, the exact relation is not known and one must resort to approximations: 
typically, good thermodynamic properties of the solid phase are obtained by using only 
the structural properties of the fluid phase as input. For instance, in the so called 
"weighted-density approximation" (WDA) [HE], one maps the free energy F[n] of the 
inhomogeneous system onto the free energy of a fluid with a smoothed density n(x), 
which is related in a non-local way to n(x). Another popular DFT is the fundamental- 
measure theory (FMT) jH] which, by enforcing dimensional crossover, provides the best 
theory invented so far for the hard-sphere solid [Zj. Recently, the FMT has been extended 
to hard-core lattice gases by Lafuente and Cuesta [S] and to soft-repulsive systems by 
Schmidt 9 . While both the WDA and the FMT work pretty well for purely repulsive 
systems, they must at present be used in tandem with perturbation theory when the 
interaction between particles has an attractive component as well JTQJ dU H21 EH Hlj- 
In this case, however, there is no guarantee that the DFT answer will be accurate. 

Recently, I have carried out (together with P. V. Giaquinta) a WDA study of 
a lattice-gas system with a realistic phase diagram (i.e., one hosting a solid, a liquid, 
and a vapour phase) jTHj, showing the effectiveness of the DFT method also for lattice 
problems. However, the theory of Ref. jT^j uses different prescriptions, in the solid 
and in the fluid phase, for the contribution of the interparticle attraction to the free 
energy, and this may represent a serious limitation if one is planning to use the same 
density functional of the bulk also for the interface problem. In particular, an interesting 
situation to analyse would be that of the solid-vapour interface at a temperature being 
only slightly below the triple-point temperature. In this case, the melting of the surface 
usually takes place, with the appearance of a thin liquid layer right at the interface (a 
successful DFT of this phenomenon for the 3D continuum case can be found in fl3j). 
As a matter of fact, we were not able to provide in J3| a sharp evidence of the surface 
melting since this would rather require that the analytic form of the (approximate) 
grand potential of the system be unique for all phases. 

In the present paper, I develop such a unique density functional for the same lattice 
model as in Ref. ^Hj by the inclusion of a convenient square-gradient term. In particular, 
the new theory leads to a genuine triple point when treating the interparticle attraction 
as a mean field, thus being a natural candidate for a practical demonstration of the 
surface- melting phenomenon in a lattice system. The toll of introducing an adjustable 
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parameter into the theory (which, as a result, partly loses contact with the original 
model) is fully compensated by the gain of a density functional which shares with the 
exact unknown one a three-minima structure in the triple-point region. It comes as a 
surprise, however, that this is not sufficient to produce an ordinary surface melting. In 
fact, an unphysical periodic phase - an artefact of the theory - is unexpectedly found 
to condense at the solid-vapour interface in place of the liquid. Only by carrying out a 
suitably constrained minimization of the density functional, aimed at washing out the 
undesired spurious phase, I finally succeed to obtain the first, at least to my knowledge, 
DFT description of complete surface melting in a lattice system. 

This paper is organised as follows. After a brief outline, in section 2, of the lattice 
DFT, I describe my system and method in section 3, while presenting results for the 
bulk and the solid surface in section 4. Finally, a summary of the main conclusions is 
given in section 5. 

2. The lattice density-functional theory: a tutorial 

The lattice DFT is a general framework for describing the statistical 

properties of systems of lattice particles under the influence of a site-dependent external 
potential e or in the presence of a self-sustained inhomogeneity (e.g. the density 
modulation of a periodic solid). The statistical description is accomplished in terms of 
the temperature (T) and chemical-potential (p) evolution of the local density n x = (c x ), 
a grand-canonical average of the occupation number c x = 0, 1 of site x in the lattice. 
n x is a lattice field which, by the Hohenberg-Kohn-Mermin (HKM) theorem, is in a 
one-to-one correspondence with the external field \x x = fi — e x , allowing one to define 
the Legendre transform F[n] of the grand potential Q with respect to the external field. 
F[n] is the inhomogeneous counterpart of the Helmholtz free energy (a more detailed 
account of the DFT method can be found in [To]). 

Given an external field fi x , one defines a sort of generalized grand potential, 



depending on two lattice fields, p and /i, which are regarded as being mutually 
independent. A minimum principle holds for saying that attains its minimum 
for a density profile n x which is precisely the one determined by fi x . Moreover, this 
minimum value is nothing but the grand potential Q for the given fi x . In practice, 
the exact profile of F[n] is not known (the ideal gas being an outstanding exception), 
and one should assign a form to F[n] which is then used for deriving, via Eq. an 
approximate grand potential for the system. 

To help the choice of F[n], one usually works with the derivatives of its excess part 
F cxc [n] = F[n] -F id [n], that is with the one- and two-point direct correlation functions 
(DCF), which, in discrete space, read as: 
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where (3 = 1/(&bT). For a fluid system with density p, spatial homogeneity imposes 
c i M = c i(p) anc ^ C i,j/N = c 2(x — y,p). Knowledge of the (two-point) DCF c 2 (x — y,p) 
allows one to obtain both the excess free energy / cxc (p) of the fluid and ci(p) (see [15J). 

Once given the fluid properties, the excess free energy of the inhomogeneous system 
can be formally expressed as an integral of the DCF through: 

f3F exc [n] = Np(3r c (p) - d(p) 5> x - p) 

X 

- J2(n x -p)(n y -p) ftX [ X aX' C V> [ nx ,] , (3) 

x,y J J 

where N is the number of lattice sites and n\ x = p + \{n x — p). At this point, 
approximations are no more eludible since the functional d^nx'] is not known. Among 
the most popular of these is the WDA |U|3], which assumes 

F exc [n] » F^ A [n] = £ nj™^) , (4) 

X 

where the weighted density n x = J2 y n y w ( x — Vi^x)- Equation (J3J) is a refinement 
of the well-known local-density approximation (LDA), -PldaW = J2x n xf exc {n x )- In 
a homogeneous system, n x is required to be constant and equal to the local density, 
implying a normalization for the weight function w; moreover, w must be such that the 
fluid DCF be also recovered in the homogeneous limit. This last condition is translated 
into a well-definite algorithm, which allows one to express n x in terms of the local 
density. For a lattice system, the details of the WDA method can be found in Ref. [T3] 
under the further simplifying hypothesis (first considered by Tarazona [4 ) of a weight 
function being a second-order polynomial in the density. 

The final step in the method is the calculation of the difference Af2[n] between the 
generalized grand potential of an inhomogeneous phase (e.g. a crystalline solid) and the 
fluid one, for equal values of T and p. One finds: 



n x 5 , 1 — 7%r 



n x \n h (1 -njln 



p 1 - p 

+ pF exc [n] - Np(3f exc (p) . (5) 

After minimizing Afi[n] over the density field, the phase-coexistence condition is 
imposed by the vanishing of this minimum difference, which implies the same pressure 
for both phases. 



3. Model and method 



As a case-study for the surface melting, I shall take a two-dimensional (2D) (rather 
than 3D) lattice-gas model with three phases. By this choice, the forthcoming analysis 
of the interface problem gets substantially simplified. Specifically, I consider the t345 
model of Ref. JH] (see Fig. 1). This is characterized by a hard-core interaction extending 
up to second-neighbor sites in the triangular lattice and a pair attraction ranging from 
third- to fifth- neighbor sites. The interaction strenghts are the same as in namely 
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v 3 = — 1.5V, t>4 = —1.2V, and t>5 = —V, with V > 0. The maximum density allowed to 
this system is p max = 0.25. The same model but with v 3 = t> 4 = v 5 = is called the 
t model. The role of the t model is that of a reference system since, in the following, 
the free energy of the inhomogeneous t345 system is going to be estimated through the 
perturbation-theory formula [H ITK] 



where the subscript refers to the t model and Av(\x—y\) = v(\x— y\) — v (\x— y\) is the 
departure of the t345 pair potential from the reference one (i.e., Av = inside the core 
while Av = v outside). In particular, it is assumed by the mean- field approximation 
(MFA) that (c x c y ) = n x n y . 

In Ref . ^5] , an approximate DCF for the t fluid is obtained by solving the Ornstein- 
Zernike relation in the mean-spherical approximation (MSA). It turns out that this MSA 
solution only exists for all fluid densities p up to p = 0.21. Beyond this threshold, the 
profile of /q xc (p) is extrapolated as a fourth-order polynomial in the density. 

The free energy of the inhomogeneous t system is expressed by the lattice 
counterpart of the WDA. The further inclusion of the mean-field attraction leads 
eventually to a density functional for the t345 model which, however, turns out to 
be insufficient to produce a liquid basin in the phase diagram [T^. In particular, the 
value of the freezing density at the critical-point temperature happens to be too small. 
To overcome this problem, one solution could be to modify the mean-field functional in 
such a way as to shift the freezing line upward in density, while leaving unaltered the 
liquid-vapour coexistence locus. To this aim, a simple square-gradient (SG) correction is 
certainly effective: it delays the appearance of the solid phase, thus indirectly promoting 
the fluid phase. In fact, the SG term has a long tradition, appearing in many popular 
phenomenological theories of the interface structure QjJ Ej . 

Indeed, a general plausibility argument for the SG correction can be provided, at 
least in case of a density field which is slowly-varying in space. As discussed in appendix 
A, the leading terms in the exact perturbative expansion of the excess free energy around 
a uniform density profile are a LDA free energy plus a SG term. However, finding out a 
direct link between this term and the microscopic-model parameters is a hard problem 
(see, however, appendix A). In fact, an even more difficult task could be to present 
a theoretical justification, within the perturbation theory, for the appearance of a SG 
correction near the mean-field attraction. In particular, including a SG correction into 
the mean-field functional might lead to overcounting some of the contributions to the 
free energy. However, if our concern is mainly to formulate a phenomenological DFT 
theory of surface melting, rather than reproducing a specific model phase diagram, the 
above objection is only of minor significance. Anyway, an ad hoc effective SG correction 
can always be regarded as just renormalizing the mean-field free energy, i.e., as a means 
for healing the crudeness of the mean-field attraction by taking into account at least 
part of the higher-order perturbative corrections. Hence, it makes sense to examine the 
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following form of free energy: 

(3F[n\ = K lnn ffi + (1 - n x ) ln(l - n x )} + Y^n x f3f^ c {n x ) 

X X 

+ 7;Y, P Av (\ x - y\) n xn y + \ ^PJM^(n x -n x+s ) 2 , (7) 

x,y * x S 

where S denotes a nearest-neighbor direction and J(p) > is, at the moment, a still 
unspecified function of the density and, possibly, also of the temperature. In the spirit 
of an entirely phenomenological approach, the virtue of (J7J) is to account for the main 
features of the t345 phase diagram, see section 4.1; then, the same theory will be asked 
to provide also a sound description of the surface-melting phenomenon. 

For the sake of simplicity, a constant value 7V is assumed for J(p). With this 
choice, the DFT free energy of the t345 system in units of k B T becomes equal, in the 
infinite-temperature, (3V — > limit, to the WDA free energy of the t system in the 
same units {i.e., the first two terms on the r.h.s. of Eq. (jZJ)). Moreover, note that the 
functional of the t model is not influenced by the SG term in Eq. (JZj). The value of 
7 follows after arbitrarily fixing the triple-point density or temperature (see the next 
section) . 

In the homogeneous limit, Eq. (JJJ) gives a generalized grand potential of 



-= plnp + (1 - p) ln(l -p)+ p(3fr(p) + \p 2 E z n (3v n - (3pp , 

n=3 



N 



(8) 



z n being the coordination number of the n-th lattice shell. The chemical-potential value 
that makes (JEJ) stationary with respect to p is: 

P 5 
(3p = ln— ci, (p) +pE z n/3v n , (9) 

1 P n=3 

where Ci i0 (p) is the one-point DCF of the t model. It has been already discussed in 
Ref. ^H] that the functional (jHJ), with (3p as in 0, shows two different minima at low 
temperature, whence two fluid phases exist, liquid and vapour, whose relative stability 
changes upon varying the density p. 

The solid phase has a triangular crystal structure that can be represented by just 
two numbers, namely the average densities and at the two sublattices A and B of 
occupied and unoccupied sites, respectively. The density functional for this phase then 
reads: 

A(3Vt^(n A ,n B ) 



N 



n A hin A + (1 - n A ) ln(l - n A ) + 3 [n B lnn B + (1 - n B ) ln(l - n B )\ 
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Upon subtracting (jSJ) from (JTUJ) and substituting Q for j3p, one finally arrives, 
after expanding the sums in (fTUj) . at the following expression for AQ: 



A(3AQ(n A ,n B ) n A 1 - n A 

= n A In h (1 - n A ) In h 3 

N p 1 — p 



, n B . l-n B 

n B m h (1 — n B ) In — 

p 1-p 



+ n A (3f™(n A ) + 3n B Pf™(n B ) - 4p/3/ exc (p) 

5 

+ 3(3v 3 n 2 A + {12(3v 4 + 6(3v 5 )n A n B + {9(3v 3 + 12/3v 4 + 6/3v 5 )n 2 B - 2p 2 £ z n (3v n 



n=3 



+ Q-yf3V(n A - n B ) 2 + ci j0 (p) - p z nfiv n (n A + 3n B - Ap) . (11) 

\ n=3 / 

The weighted densities n A and n B are evaluated from the values of n A and n B with 
the formulae appearing in appendix A of Ref. ^H]- In practice, the calculation of n A 
and n B is carried out in parallel with that of n A and n B , i.e., by the same iterative 
procedure which determines the minimum of AQ. This could be a steep est- descent 
dynamics or, equivalently, the self-consistent solution of the equations expressing the 
stationarity condition for AQ. These equations are easily found to be: 

n~ A l = l + p - exp \ Cl (p)-pJ2 Zn(3v n + /?/ exc (n A ) + n A (3r c \n A )— A - 



n=3 

f cxc / 



- 3n B (3f exc '(n B )^ + Q(3v 3 n A + (12(3v 4 + Q(3v 5 )n B + 12^V(n A - n B ) 
dn A 



n 



B 



1 + — exp ( Cl (p) - p £ Znf3v n + Pr c (n B ) + \n A (3r c \n A )^ 
P I ^3 3 9n B 

+ n B (3r c \n B )^ + (Af3v 4 + 2(3v b )n A + (Qf3v 3 + 8f3v 4 + A(3v 5 )n B + A 1 (3V{n B 

OUR 



I refer the reader again to ^3] for more information about the technicalities of the 
minimization procedure for AQ. 

4. DFT results 

Taken the functional (jllj) for granted, I now review the results for the bulk properties 
of the system. Next, I refer on the attempt of using the same functional also for the 
solid-vapour interface. 

4-1. The bulk 

First, I summarize the results of the WDA theory for the t model. The coexisting fluid 
and solid densities are found to be pf = 0.1335 and p s = 0.1686, respectively [T^]. The 
actual values, obtained through the grand-canonical Monte Carlo (MC) method ^H] , are 
instead pf = 0.172(1) and p s = 0.188(1), indicating that the thermodynamic stability of 
the t solid is overestimated by the theory. 
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Moving to the attractive model, I first seek for the minimum of JSJ), with f3fi as 
in Eq. Q. For low enough values of the (reduced) temperature t = k B T/V, there 
are in fact two distinct points of minimum, the deeper one being associated with the 
thermodynamically stable phase. The locus of points (p, t) where the two minima have 
equal depth is the liquid- vapour coexistence line (the crosses in Fig. 2). In particular, 
the critical point falls at p cr = 0.068(1) and t CT = 0.778(1). 

As far as the solid-fluid equilibrium is concerned, I first adjust the quantity 7 so as 
to obtain a triple point at a selected temperature t tv - It turns out that, in order for t tr 
to be e.g. exactly 0.7, the value of 7 in (jllj) has to be 0.25191(1). A bit larger value, 
i.e., 7 = 0.28214(1), shifts the triple-point temperature down to 0.6. 

The complete DFT phase diagram of the t345 model is plotted in Fig. 2 for three 
values of 7, namely (MFA), 0.25191 (t tr = 0.7), and 0.28214 (t tr = 0.6). The MC 
simulation data of Ref. ^HJ are shown as asterisks. It clearly appears that the functional 
(I11J) yields a realistic phase diagram only in a small window of 7 values. In Ref. |15j . 
the predicted DFT phase diagram was equally good (Fig. 2 of [E], open circles), but we 
were obliged to use distinct DFT functionals for the solid and the fluid. A minor defect 
of the present theory is the unphysically large (> 0.25) value of the melting density at 
low temperature. In fact, though and % are constrained to be smaller than 1, their 
combination n s = (n^ + 372B)/4 is not under control during the functional minimization. 
In any event, the agreement of the DFT with MC remains mainly qualitative, also 
because of the not excellent quality of the WDA theory for the reference t system. 

As already mentioned, the major advantage of (fTTj) over the functional considered in 
Ref. [T3] lies in the fact that it has a unique expression for all phases of the system. This 
is evidenced e.g. by its profile at the triple-point temperature: for 7 = 0.25191 (i.e., 
t tr = 0.7), the three-minima structure of Afl(riA_, Ub) is shown in Fig. 3, as projected 
onto the % = plane in the 3D space spanned by n^n^, and AQ. Upon moving 
away from the triple point, the relative depth of the minima changes, thus making a 
particular phase of the system more stable than the others. 

4-2. The solid-vapour interface 

A density functional with three minima, like that at Eq. (jllj) . gives the rather unique 
opportunity to describe on an equal footing, i.e., within the same theoretical framework, 
both the bulk and the surface statistical properties of a many-particle system. Among 
these, a prominent place is certainly held by the surface-melting phenomenon, which 
opens a sort of imaginary window on the structure of the underlying generalized 
grand potential: roughly speaking, it is the liquid state claiming visibility before 
thermodynamics gives its allowance. 

Let me consider, for instance, a linear interface running along a nearest-neighbor 
direction, called X. This interface breaks the translational symmetry along the 
perpendicular, vertical direction Y, thus causing the sublattice densities to get a Y 
dependence. Very far from the interface, the densities recover the bulk values, being 
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those of the solid for e.g. Y <C and those of the coexisting vapour for Y ^> 0. The 
horizontal layers are labelled with an index A, which increases upon going from solid to 
vapour, being zero at the "centre" of the interface. To be specific, those layers where 
particles are preferentially hosted in the solid have an odd A value. At variance with the 
bulk case, three sublattices are to be distinguished now, since different density values are 
generally expected at the even and at the odd interstitial sites. I call C the sublattice 
formed by the interstitial sites pertaining to the odd layers, and B the other. Finally, 
A is the triangular sublattice that is occupied in the T = solid. 

Then, a rather straightforward adaptation of (fTUj) to an interface geometry yields 
eventually a surplus £[n] = 2f3/S.Q\n]/Nx of grand potential per surface particle, due to 
the interface, equal to: 



E 

A odd 



1 UA ' X , ft M 1 ~ n A,\ , , n C,X . n M 1 ~ n °> X 

n A ,A In h (1 - n A ,x) In — h n c ,\ In h (1 - n c ,\) In - 
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+ C l,o(p) - pJ2 Z nP l 
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E ( n Ax + n c,x - 2p) + 2 ( n B,x ~ p) 



.X odd 



A even 



+ E [n A ,xf3fr^A,x) +n c ,x(3fr(nc,x) -2p(3fr(p)} 

X odd 

+ 2 £ K A /?/ exc (n BiA )-p/3/ cxc (p)] 

A even 
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A odd 
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A even 



n A,A E n vPAv(\x - y\) +n c ,x E ^j//5Ai;(|a; - y\) - 2p 2 ^ z n 

y\x<=A,X y\x£C,X n=3 

5 

^b,a E n y /3Av(\x - y\) - p 2 ^ z n (3v rt 

y\x€B,X 



n=3 



1 A odd 



E ( n AA - n x+5 ) 2 + ^ (nc,x - n x+& ) 

&\x&A,X S\xeC,X 



+ 2 E E i n B,x - n x+s f 

X even (5|x€B,A 

To save space, some of the sums in the above Eq. ()1HJ) are left as indicated; however, 
expanding them is not particularly difficult, it is just a matter of carefully looking at 
the lattice geometry. 

I have considered the case 7 = 0.25191 and a number of temperature values in the 
range from 0.6 to t tr = 0.7. A slab of interface generally consisted of 160 layers, from 
A = —80 to A = 79, with fixed values of the sublattice densities outside this range. The 
extension of the slab in the X direction is virtually infinite. The interface is prepared 
by always assuming an inverted-exponential, (1 +exp[(A + 1/2)//]) -1 modulation [TH] 
which is then relaxed until a minimum of E[n] is found. Disappointingly, however, the 
final looking of the interface was never as expected. A typical outcome is depicted in 
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Fig. 4, which refers to t = 0.65. We can see that the mean density of a layer, ri\, shows 
a curious ten- layer-long periodic motif for A > 0, which interrupts at the right border 
of the picture for the existence of a fixed boundary. Even more serious is the fact that 
the minimum of £[n] turns out to be negative, signalling that there is something wrong 
with this interface. 

To clarify such things better, I have cut out a piece of this periodic motif and made 
it infinite in both Y directions. The resulting phase does, in fact, overcome in stability 
both the solid and the vapour. Obviously, such a weird evidence is totally unphysical, 
being a spurious result of the DFT, as founded on Eq. (JT3J). It is worth noticing that, 
in this periodic phase, the local values of the system density are different at the three 
sublattices A, B, and C, i.e., this phase can only exist in an interface. Nonetheless, 
the very occurrence of this oddity inevitably casts a shadow on the real quality of the 
functional (jJJ), which should then be re-examined in spite of its capability of accounting 
for the bulk properties of the t345 model. 

Two different routes are now opening to our consideration: the first is simply 
to reject the functional (J7|) as unphysical. Another solution is however viable, which 
is to insist on (jJJ), trying to see whether it is possible to unveil the normal surface- 
melting behaviour by properly obstructing the manifestation of the undesired phase. 
Obviously, this is not completely orthodox since, in a sense, a sort of external field is 
being introduced into the problem. 

If we decide to pursue this second route anyway, one can try the following. Rather 
than seeking for the absolute minimum of S[n], the search is restricted to just those 
fields n x that, at a particular C site, take the same value as in one of the two closest 
B sites. Precisely, I require that nc,x = tib,x-i, a choice which is consistent with the 
boundary values but incompatible with the unphysical periodic phase. Whether this 
is sufficient to promote the appearance of the liquid at the middle of the interface is a 
matter of debate that can only be fixed numerically. 

In order to carry out the constrained minimization of E[n] along the guidelines 
presented above, Eq. ([T3J) must be modified by simply substituting tib,\-\ for na,x 
everywhere it occurs. A more subtle, but absolutely crucial, technical point is relative 
to the weighted densities: it follows from their very definition that there is actually no 
simple relation between nc,x and n^-i, even when ric,x = Ub,\-i- This is due to the 
fact that, owing to the Y dependence of the sublattice densities, the neighborhood of a 
C site does not look the same as that of a B site, and this causes nc,x and ub,\-i to be, 
at least in principle, different. In the end, the amended E[n] functional reads as: 
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n=3 / LA odd A even 

+ £ [n Ax f3fT(nA,x) +n B ,x-iPfr(nc,x) -2p(3fT(p)} 

X odd 



Lattice density-functional theory of surface melting 



11 



+ 2 J2 [n B ,x(3fr(nB,x)-p(3fr(p)} 

X even 

+ o H W.A [2/^ 3 (n AA -2 + 7U,a + ^4,a+2) + 2/3u 4 (2n Bj A_3 + ri B)A -i + 2n BjA+ i + n B ,A+3) 

^ A odd 

+ 2(3v 5 {n BtX -3 + n B ,\-i + n B<x+3 )} ~ p 2 J2 z n(3v n \ 

n=3 J 

A even 

+ (3v A (n B) x-A + n AtX -3 + 4n B ,x-2 + 2?2a,a-i + 2n BiA + ru,A+i + 4n BiA+2 + 2ru,A+3 + n BtX +4) 
+ (3v 5 (n B ,\-4 + n A ,A-3 + ^b,a-2 + 2n BiA + n AtX +i + n B ^+2 + ^a,a+3 + n Bt \ +4 )] 

~ \p 2 ^2 z n(3v n \+'y(3v\ [3(n A ,x ~ ^b,a-i) 2 + {n A ,\ ~ n B ,x+i) 2 

Z n=3 J KXodd 

+ [i n B,x - n AX -i) 2 + 2(n BjX - n B ^-2) 2 + (n B ,x - n AX +i) 2 > ■ (14) 

A even ) 

I omit to indicate the explicit expression of the weighted densities and of their 
derivatives, since they would need too much space to be specified here. 

Using the functional (JHJ), I have carried out a new series of minimizations, for 
various temperature values. With great pleasure, the thermal evolution of the interface 
structure is now as expected in a 3D system with a complete surface melting, see Fig. 5, 
though dissimilar to the behaviour of the actual (2D) t345 model where the layering of 
the liquid is much more pronounced, cf. Fig.9 (top) of Ref. [T3]. Moreover, the present 
evidence is also superior to that found in jTH| by using different functionals for the fluid 
and the solid. There, the thickness of the molten layer was too small even very close to 
the triple-point temperature. 

The slab used in the present optimizations was a hundred layers wide, from A = —50 

to A = 49, but only a slice of this is shown in Fig. 5. In Fig. 6, I report the detailed 

profile across the interface of the sublattice densities for two distinct values of t close 

to t tT , namely t = 0.69 and t = 0.699. By looking at Fig. 5, we see that the width 

of the liquid layer grows with regularity as the triple-point temperature is approached 

from below, corresponding to the liquid phase becoming more and more stable. To 

estimate this width as a function of temperature, I use the following three-parameter, 

double-exponential fit: 

n A - p\ pi- p 

n A \=p-\ — ; — — - H — ; — — , for odd A only ; 

' H l + exp[(A + l/2 + L)/Zi] l + exp[(A + l/2-L)// 2 ] ' 

kb- Pi pi-p , . . 

n B \ = p-\ p— ; — — H — ; — — , for even A only ; 

' l+exp[(A + l/2 + L)/Zi] 1 + exp [(A + 1/2 - L)/Z 2 ] ' *' 

n>c,\ = n Bt \-i, for odd A only. (15) 

In Eq. ()15p. p\ is the density of the (metastable) liquid coexisting with the vapour at the 
given temperature t < t tT , 2L is the nominal width of the liquified part of the interface, 
while l\ and Z 2 are a measure of the width of the two transition regions being centred at 
— L — 1/2 and L — 1/2, respectively. The values of the fitting parameters are obtained 
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by substituting the ansatz (115)) into Eq. (|14j) and requiring it to be minimum. It turns 
out that this minimum value is never appreciably far from the absolute minimum of 
S[n] as being determined independently. 

For t = 0.65, 0.69, and 0.699, 1 find L = 2.67, 4.61, and 8.60, respectively (cf. Fig. 5). 
Hence, while at low temperature it is more convenient for the solid-liquid and liquid- 
vapour interfaces to stay bound together, slightly below the triple point the free-energy 
cost of two well-separated interfaces is smaller. In some phenomenological theories of the 
surface melting, this effect is mimicked through the introduction of an effective repulsion 
between the two interfaces. Furthermore, while l\ ~ 1.0 shows only a small dependence 
on temperature, the increase of I2 with t is more conspicuous, being 1.42, 1.81, and 2.86 
for the above cited t values. At these same temperatures, the minimum E[n] turned out 
to be 1.10162,0.74472, and 0.65504, respectively. In particular, the latter value would 
practically account for the overall cost of two separate solid-liquid and liquid-vapour 
interfaces at t = t tT . I do not try to deduce from the above numbers an empirical 
relation between 2L and t, e.g. with the aim at resolving a logarithmic from a power- 
law behaviour. In fact, owing to a necessarily imperfect estimate of the coexistence 
conditions within the DFT, it was actually impossible to approach the triple point more 
closely than 0.001 in t. This notwithstanding, I do not find any reason to doubt that L 
would grow to infinity for t — > t^ T . 

Summing up, the surface melting is observed, within the functional (fTTj) . only at 
the price of artificially restricting the symmetries that are available to the density field. 
Otherwise, in fact, a spurious periodic phase will wet the solid, which is actually more 
stable than the solid and vapour themselves. 
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5. Conclusions 

The usual realm of applications of the classical DFT comprises hard-core systems 
and soft-repulsive ones. Besides those cases, the only existing general DFT method 
is, at present, the density- functional perturbation theory which, however, can be 
quantitatively inaccurate or, even, predict erroneous phase diagrams. 

In this paper, I have used the lattice analogue of the perturbative DFT in order to 
analyse the phase behaviour and the interface structure of a 2D lattice-gas model with a 
pair interaction consisting of a hard core plus an attractive tail. From previous studies, 
this model is known to exhibit a solid, a liquid, and a vapour phase. 

In order to obtain a phase diagram with three phases, I have considered a 
phenomenological density functional which, besides treating the interparticle attraction 
as a mean field, also contains a tunable square-gradient correction. In fact, the 
very unique request to this theory was to yield a generalized grand potential for an 
inhomogeneous phase that could exactly turn, in the homogeneous limit, into the fluid 
functional. This requirement is crucial for getting a consistent description of the interface 
between the solid and a fluid phase. 

After producing a good theory for the bulk, I have adapted the same density 
functional to the geometry of a solid-vapour interface. The aim was to observe, under 
equilibrium conditions, the melting of the solid surface, i.e., the growth of a liquid layer 
between the solid and the vapour as the triple point is approached from low temperature. 

Surprisingly, the minimization of the interface functional unveils instead a periodic 
phase which is totally unphysical, being actually an artefact of the theory. In order 
to recover a more conventional behaviour, I have carried out a constrained functional 
minimization which, while not affecting the liquid, makes it unlikely for more exotic 
phases to appear. By this stratagem, and nothwistanding the partial arbitrariness that 
is implicit in such a conditioned minimization, in the end I arrive at a fairly realistic 
description of the complete surface melting of a lattice system - a result which, in my 
opinion, is an interesting illustrative example of the risks which may accompany the 
use of a phenomenological DFT functional. Anyway, I do not think that the ad hoc 
expedient that was considered here can be of help for investigating, in the absence of a 
good-quality functional, more delicate features like e.g. layering phase transitions and 
roughening. 
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Appendix A. Estimate of the square-gradient correction 

In this appendix, I show how to adapt to a discrete space an argument originally due 
to Evans about the estimate of the "optimal" SG correction to a LDA functional pQ. 
Let the following excess free energy be considered, 

F^[n]=Y j nJ exc (n x ) + Wj{n x )Y J (n x -n x+5 )\ (A.l) 

x x S 

which is a LDA free energy supplemented with a SG term. This form of F exc [n] is not 
precisely the one which I work with in section 3 for demonstrating the surface-melting 
phenomenon. Actually, my intention here is not to provide a full justification of Eq. ([7)1; 
rather, I just want to give a flavour of the physical status of an SG correction. 

To state the problem precisely, the function J(p) in Eq. ((7J) will be calculated from 
the requirement that at least the second-order expansion of (jA.l|) around a fluid of 
density p may agree with the exact one, 

F cxc [n] = N P r c {p) - I Cl (p) £ An x - ±- £ c 2 (x - y, p)An x An y + ... (A.2) 

P x r x,y 

The field An x = n x — p is a measure of the (admittedly small) deviation of the system 
from homogeneity. 
Since 

c 1 (p) = -pr c (p)- P [3r c '(p) J (a.3) 

the expansion of the LDA term in (|A.1J) is as follows: 

E^f^K) = N P r c (p) - i Cl (p) £ An, - ^(P)i: Anl + ... (A.4) 

To proceed further, it is convenient to work in Fourier space. I take the Fourier transform 
of a lattice field f x to be f q = J2x fx exp(— iq-x) (see the precise definition of the q vectors 
in Ref. 17 of [15J). The convolution theorem first yields: 

E A ™1 = I E ^An- q ; (A.5) 

x Q 

similarly, the fully quadratic, leading SG correction becomes as follows: 

\j{p) EK - n x+s ) 2 = J{p) £ W{x - y)An x An y = J{p)^- £ W q An q An_ q . (A.6) 

1 x,S x,y iV q 

The function W in Eq. (|A.6|) is lattice-dependent; for e.g. a d- dimensional hypercubic 
lattice with unit lattice constant, W(x — y) — 2d5 x>y — and W q = 2d[l — 

(1/cOELiCosgJ. Finally, 

E c ^ x ~ V> p)^n x An y = 4 E c ~2(?, p)An q An- q . (A.7) 

x,y iV q 

Upon comparing Eqs. (|A.4jl . (|A.5|) . and (|A.6|) with Eq. (|A.7jl . one would be led to 
identifying c[(p) — 2f3J(p)W q with c 2 (q,p). However, it is clear that there is no chance 
for these two quantities to be exactly equal, since J(p) does not depend on q. To fix this 
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problem, one should make some further assumption on the An x . For instance, one can 
require that, besides being small, An x also varies slowly in space. In this case, only the 
long- wavelength Fourier components of An x are important, implying that what really 
matters in fact is the matching of the above two quantities in the small-g limit. In 
particular, the LDA free energy (i.e., J(p) = 0) does exactly reproduce only the q = 
Fourier component in Eq. (|A.7|L since c[(p) = c 2 (0, p). 
At second order in q, 

p) = c 2( x > P) ( x ~ ' l Q ' x ~ \(l ■ x f + • • •) • ( A - 8 ) 

The first term J2 X c i(x, p) = 62(0, p). The second term is zero for a lattice with inversion 
symmetry (i.e., one with C2(x,p) = C2(—x,p)). The third term is: 

~^Y c 2(x,p)(q-x) 2 = -^q 2 Y,x 2 c 2 (x,p)(q-x) 2 = - ^q 2 ^ x 2 c 2 (x , p) = d 2 (p)q 2 , (A.9) 

having denoted unit vectors by a hat. Observe that, in the last step, the axial symmetry 
of the lattice with respect to each coordinate axis has been assumed. Hence, the 
expansion of c 2 (g,p) in powers of q begins as c 2 (0,p) + d 2 (p)q 2 . 

Since W q = q 2 + o(q 2 ) for a hypercubic lattice, the expression of J(p) for this lattice 

is: 

Ap) = -jpMp) = J^E^2(^,P). (A.10) 

For a triangular lattice, W q = (3/2)q 2 + o(q 2 ) while o^G ) is like in Eq. (|A.9j) with d = 2. 
Hence, J(p) is still given by Eq. (jA.10|) . but with d = 3 (the rationale being that, for 
the triangular lattice, the nearest-neighbor sites are six, like for the cubic lattice). 

To conclude, I make some comments on the result ()A.10|) . For the t model, the 
non-zero, core values of C2(x,p) are all negative, hence J(p) < 0. This means that the 
SG correction is meaningless for this model. Different is the case of the t345 model, 
where the positive tail of the DCF would eventually make J(p) > 0. However, the exact 
C2(x,p) is not known, and this actually makes Eq. (jA.10|) rather useless. 
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Figure Captions 

Fig. 1 : The triangular-lattice model under consideration (schematic): excluded sites 
(x) and attractive sites (numbered dots) are separately shown for a particle sitted 
at the centre of the picture (large black dot). The attractive interaction reaches 
the fifth neighbours, no interaction is felt beyond this distance (small dots). The 
Hamiltonian is of the form H = J2i<j U (N — j\) c i c j> where v(\i — j\) = +00 when 
simultaneous occupation of sites i and j is forbidden. In the text, v n denotes the 
value of v(\i — j\) for a pair of n-th neighbours, whereas z n is the number of n-th 
neighbours (z% = Z2 = z 3 = 6, Z4 = 12, and Z5 = 6). In this work, the interaction 
strengths were v 3 = —1.5V, V4 = —1.2V, and v 5 = —V (with V > 0). For these v n , 
the phase diagram of the lattice model is the canonical one, with a solid, a liquid, 
and a vapour phase (see Ref. [IE]). 

Fig. 2 : Phase diagram of the t345 model as drawn from the density functional (fTTj) . 
using the t model (MSA+WDA) as a reference and representing the interparticle 
attraction as a mean-field perturbation. Various solid-fluid coexistence lines are 
shown: 7 = (A and dotted lines) [T3], 7 = 0.25191 (O an d continuous lines), 
and 7 = 0.28214 (□ and continuous lines). The liquid- vapour coexistence locus is 
marked with crosses. For comparison, I have reported as asterisks joined by dashed 
lines some MC data points for the 48 x 48 lattice (the errors affecting these points 
are of the same size as the symbols). The arrows pointing downwards mark the 
WDA densities of the coexisting fluid and solid in the t model. The other arrows 
point to the MC estimates for the same quantities. 

Fig. 3 : The function AQ(n^, Ub) at Eq. (jllj) is plotted vs. ha, for various n% values 
ranging from to 1 (both sublattice densities are incremented in steps of 5 x 10 -3 ). 
Here, 7 = 0.25191 and t = t tr = 0.7. Straight lines are drawn through the 
points as a guide to the eye. The vapour and liquid densities are p v = 0.0248 and 
Pi = 0.1282, respectively. The coordinates of the "solid" minimum are = 0.94767 
and hb = 0.02765, corresponding to a solid density of p s = 0.2577. 

Fig. 4 : DFT profile of the mean density n x at the A-th layer, defined as (nA,\ + n c ,\)/2 
for odd A; as ns,\ for even A. What is actually plotted is the outcome of the 
minimization of (|T3*j) for 7 = 0.25191 and t = 0.65. During the minimization 
process, it is noted that the vapour disappears progressively in favour of a novel 
phase with a period of ten layers. The automatic search of the minimum AQ comes 
to a stop when this periodic phase reaches the border of the slab. 

Fig. 5 : DFT results from the minimization of functional (fT4*)) for 7 = 0.25191 (t tT = 
0.7): the optimum mean-density profile across the solid- vapour interface is plotted 
for three values of t, namely 0.65, 0.69, and 0.699. To help the eye, in each panel 
a continuous line has been drawn through the points. As is quite clear from the 
pictures, the typical surface-melting behaviour of a 3D simple fluid is derived from 

Eq.dU. 
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Fig. 6 : Layer- by-layer evolution of Ua,\ (for odd A only), Ub,\ (for even A only), and 
n c,x — n B,\-i (f° r odd A only) for the two same interface profiles of highest 
temperature as in Fig. 5. Near A = 0, the plateau of ua,x and the maximum of 
both ub,\ and nc,\ are clear signs of a local liquid-like behaviour. 
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Fig. 1 (S. Prestipino, Lattice DFT...) 




Fig. 2 (S. Prestipino, Lattice DFT...) 
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Fig. 4 (S. Prestipino, Lattice DFT...) 
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Fig. 5 (S. Prestipino, Lattice DFT...) 
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Fig. 6 (S. Prestipino, Lattice DFT ...) 



